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I.  INTRODUCTION 


A  protuberance  in  a  boundary  layer  whose  heiqht  is  of  the  order  of  the 
boundary  layer  thickness  would  be  expected  to  cause  large  changes  in  the 

neighboring  flow  field.  Experimentally  it  has  been  found  that  significant 
perturbations  in  the  flow  are  also  induced  downstream  of  the  protuberance  for 
hundreds  of  protuberance  heights.  As  determined  by  flow  visualization  and 

supported  by  hot  wire  measurements,  these  perturbations  take  the  form  of 
streamwise  vortices  or  concentrations  of  streamwise  vorticity.  A  theoretical 
model  for  this  persistent  effect  does  not  exist;  in  this  report  numerical 
modelling  is  employed. 

There  are  other  flow  fields  where  the  persistence  of  streamwise  vorticity 
has  been  found.  In  wind  tunnel  experiments  on  flat  plates,  non-uniformities 
in  screens  far  upstream  of  the  plate  have  caused  significant  spanwise  velocity 
components.  Small  non-uniformities  on  the  leading  edge  of  plates  or  on  steps 
of  bodies  of  revolution  also  induce  persistent  streamwise  vortices.  Such 
vortices  are  also  formed  in  the  boundary  layer  on  the  leeward  side  of  a  body 
of  revolution  at  moderate  angle  of  attack  and  can  persist  to  the  base;  as 
angle  of  attack  increases,  the  vortices  are  shed  from  the  body.  These  vorti¬ 
ces  can  influence  the  Magnus  effect  on  a  projectile. 

The  persistence  of  streamwise  vortices  is  found  for  both  laminar  and 

turbulent  boundary  layers.  For  the  latter,  spanwise  non-uniformity  in  the 
transition  process  can  induce  them  in  addition  to  the  mechanisms  mentioned 

above.  Streamwise  vorticity  may  also  be  involved  in  large-scale  structures  in 
turbulent  boundary  layers.  A  matter  of  practical  concern  when  streamwise 
vortices  exist  in  turbulent  boundary  layers  is  that  they  increase  the  heat 
transfer. 

The  above  ideas  are  discussed  by  Morkovin1,  who  also  emphasizes  the 
unsatisfactory  state  of  the  observations  and  theory  for  streamwise  vortices. 
The  main  conclusions  of  the  present  work  are  also  discussed. 

To  study  the  effect,  the  simplest  situation  is  chosen  in  which  streamwise 
vortices  are  found:  flow  behind  a  protuberance  imbedded  in  a  low  speed  lami¬ 
nar  boundary  layer  on  a  flat  plate.  Also,  more  experimental  data  is  available 
for  this  case  than  any  other;  even  so  it  is  not  sufficient  for  our  purposes. 
The  objective  of  this  work  is  to  determine  if  numerical  simulation  of  the  flow 
field  will  show  the  persistence  and  other  features  of  the  flow  field  that  are 
observed  experimentally. 

A  schematic  of  the  flow  field  is  shown  in  Figure  1.  Qualitatively,  the 
same  flow  field  features  are  found  over  a  wide  range  of  conditions,  regardless 
of  whether  the  boundary  layer  is  laminar  or  turbulent,  or  whether  the  external 


1.  Morkovin,  M. V.}  "Observations  on  Streamvise  Vortices  in  Laminar  and 
Turbulent  Boundary  Layers ",  NASA  Contractor  Report  159061 ,  April  1979. 
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flow  is  high  speed  or  low  speed.  Sedney2  has  surveyed  the  effects  of  small 
protuberances  on  boundary  layer  flows;  for  additional  details  of  the  flow 
field  and  references  describing  them  see  Reference  2.  As  sketched  in  Figure 
1,  the  boundary  layer  separates  upstream  o-'  the  protuberance,  forming  a  horse¬ 
shoe  vortex  which  stretches  around  the  obstacle.  Spiral  vortices  rise  up  from 
the  surface  of  the  plate  in  the  near  wake,  forminq  twin-vortex  filaments  which 
trail  downstream.  Smoke  visualization  studies  in  'ow-speed  boundary  layers 
show _  that  the  streamwise  vortices  produced  by  the  interaction  process  can 
remain  steady  and  distinct  for  hundreds  of  protuberance  heights  downstream  of 
the  obstacle.  The  persistance  of  these  vortices  is  quite  remarkable,  even 
more  so  when  it  occurs  in  turbulent  boundary  layers,  considering  the  high 
diffusivity. 

Two  numerical  approaches  for  modeling  the  downstream  flow  field  are 
employed  in  the  present  study.  The  first  uses  3-D  boundary-layer  approxima¬ 
tions;  the  second  makes  use  of  a  3-D  boundary-region  approximation  which 
includes  the  viscous  crossflow  (diffusion)  terms  neglected  in  the  boundary- 
layer  theory.  The  measurements  of  Tani,  et  al.3’  taken  ten  protuberance 
heights  downstream  of  the  obstacle,  are  used  to  construct  the  initial  data  for 
the  calculations.  Quantitative  agreement  between  the  numerical  simulation  and 
measurements  will  depend  on  the  accuracy  with  which  we  can  fit  the  initial 
data,  which  is  obtained  from  experimental  measurements  and  properties  of 
boundary  layers.  The  uncertainties  in  the  fitting  wi'l  be  described. 

The  major  conclusion  is  that,  using  the  boundary-region  approximation,  we 
are  able  to  reproduce  the  qualitative  features  of  the  flow  field  far  down¬ 
stream  of  the  obstacle,  including  the  persistence  of  streamwise  vorticity;  but 
we  do  not  succeed  using  the  3-D  boundary-layer  approximation. 


II.  GOVERNING  EQUATIONS 

For  some  range  of  the  parameters  involved,  si. ch  as  Reynolds  number, 
protuberance  dimensions  and  geometry,  the  flow  is  sensibly  steady;  see  Refer¬ 
ence  2.  Therefore  we  use  the  steady,  3-D,  incompressible  Navier- Stokes  equa¬ 
tions.  * 


u  u 


x 


+  v  u 

y 


+  W  U2 


"<uxx  + 


yy 


uzz>’ 


(1) 


2.  Sedney,  R.,  "A  Survey  of  the  Effects  of  Smalt  Protuberances  on  Boundary- 
Layer  Flows",  AIAA  Jowmal,  Vol.  11,  No.  6,  June  1973,  pp.  782-792. 

3.  Tani,  I.,  Komoda,  H. ,  and  Komatsu,  Y.,  "Boundary-Layer  Transition  by 
Isolated  Roughness",  Aeronautical  Research  Institute  Report  No.  375, 
University  of  Tokyo,  November  1962. 

*Definitions  of  symbols  are  given  in  the  LIST  OF  SYMBOLS  section. 
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(4) 


w  are 


the  density,  p  the 
dimensional  velocity 


where  subscripts 

pressure,  p  the  absolute  viscosity;  and  u,  v, 

components  in  the  directions  of  the  streamwise  coordinate,  x,  the  spanwise 
coordinate,  y,  and  the  normal  to  the  plate,  z,  respectively  (see  Figure  1  for 
coordinate  system).  Equations  (l)-(4)  are  sufficient  for  describing  the 
complete  protuberance  flow  field,  but  would  be  very  difficult  to  apply  in 
practice.  If  we  concentrate,  instead,  on  the  downstream  flow  field,  then  we 
suspect  that  if  we  initialize  the  calculations  "far  enough"  downstream,  we  can 
use  the  boundary- layer  or  boundary-region  equations  to  describe  the  flow 
field. 

The  3-D  boundary  layer  equations  can  be  derived  from  Equations  (l)-(4) 
through  either  an  order-of-magnitude  analysis,  or  by  formally  applying  limit 
processes  at  high  Reynolds  numbers.  The  resulting  equations,  in  non-dimen¬ 
sional  form,  are 


U  U, 


V  Uy  + 


w  uz  = 


-Pv  +  U 


11 


(5) 


U  Vv  +  V  Vw  + 


w  Vz  -  -PY 


+  V 
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(6) 


Ux  +  Vy  +  Wz  =  0 

where  the  non-dimensional  variables  are  defined  by 

X  =  x/L,  Y  =  y/L,  U  =  u/u  ,  V  =  v/u 


Z  =  z  Re1/2/L,  W  =  w  Re172^,  P  =  p/(pu J) 


(7) 

(8) 
(9) 


with 


Re  =  puJ-/p. 


(10) 


The  freestream  velocity  is  u^  and  L  represents  a  reference  length,  taken  to  be 
1  cm  in  the  present  problem. 

The  3-D  boundary-region  equations,  on  the  other  hand,  cannot  be  formally 
derived  from  Equations  (l)-(4),  to  the  best  of  our  knowledge.  Rather,  they 
represent  an  ad  hoc  approximation  that  has  been  used  in  some  other 
investigations  to  treat  problems  where  crossflow  diffusion  effects  were  not 
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negligible.  They  were  first  used  by  Kemp4  for  the  flow  over  two  intersecting 
flat  plates,  the  corner  problem,  and  have  been  employed  extensively  by  many 
authors  for  that  problem.  The  boundary- region  equations  can  be  derived,  using 
familiar  limiting  processes,  for  the  corner  problem,  but  not  for  the  present 
problem  or  many  others  to  which  they  have  been  appliec.  Therefore,  to  account 
for  the  expected  crossflow  diffusion  effects,  we  adopt  the  boundary-reqion 
approximation  expressed  by 


U  Ux  +  V  Uy  ♦  W  Uz  =  -Px  ♦  Uzz  ♦  Re'1  Uyy 

(11) 

t  1  V,  t  IUZ  =  -P,  .  Vzz  +  Re'1  Vyy 

(12) 

Ux  +  VY  +  Wz  =  0, 

(13) 

where  the  unusual  form  of  the  diffusion  terms  on  the  right  hand 
from  the  way  in  which  z  is  non-dimensional ized. 

side  results 

We  can  write  a  single  set  of  equations  that  expresses  both 
layer  and  boundary-region  approximations 

the  boundary- 

U  Ux  +  V  Uy  +  W  Uz  =  -Px  -  e  Re”l  Uyy  +  Uzz 

(14) 

U  Vx  +  V  Vy  +  W  Vz  =  -Py  -  e  Re-1  VyY  +  Vzz 

(15) 

Ux  +  Vy  +  Wz  =  0, 

(16) 

where  e  =  0  gives  the  3-D  boundary- 1  ayer  equations  and  e  =  1  gives  the  3-D 
boundary-region  equations.  The  boundary-region  approximation  includes  one 
additional  crossflow  diffusion  term  in  each  momentum  equation  that  is  neglect¬ 
ed  in  the  3-D  boundary- layer  theory.  Both  approximations  neglect  the  normal 
momentum  equation,  assuming  that  normal  pressure  gradient  Pz  =  0;  model  equa¬ 
tions  with  Pz  *  0  could  also  be  developed.  The  boundary  conditions  applicable 

to  Equations  (14 ) - ( 16 )  will  be  described  in  Section  V.  When  the  external  flow 
is  uniform,  as  is  assumed  here,  P^  =  Py  =  0. 


III.  TANI  TEST  PROBLEM 

We  have  attempted  to  apply  Equations  (i4)-(16)  to  describe  the  flow  field 
downstream  of  a  typical  small  cylindrical  protuberance.  Experimental  hot  wire 
anemometry  measurements  taken  by  Tani ,  et  al.3  establish  the  presence  and 
persistence  of  streamwise  vortices  in  the  downstream  flow  field  for  various 


4.  Kemp. N .  "The  Laminar  Three-Dimensional  Boundary  Layer  and  a  Study  of  the 
Flow  Past  a  Side  Edge",  M.Ae.S.  Thesis,  Cornell  University,  1951. 
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free  stream  conditions.  The  test  problem  treated  here  corresponds  to  Tani's 
experimental  arrangement  with  a  cylindrical  protuberance  (k  =  0.25cm),  where  k 
is  the  protuberance  height,  mounted  normal  to  a  flat  plate  placed  in  a  wind 
tunnel;  see  Figure  1.  The  protuberance  height  is  equal  to  its  diameter  and 
the  freestream  velocity  is  520cm/s,  with  a  free  stream  turbulence  level  of 
approximately  0.2%.  The  protuberance  leading  edge  is  located  at  xk  =  60  cm. 

For  these  conditions  the  protuberance  height  is  approximately  39%  of  the  local 
boundary- layer  thickness,  6,  at  x  =  xk,  based  on  the  99%-velocity  thickness. 

The  Reynolds  number  based  on  xk,  with  v  =  u/p  =  0.15cm2/s  is 


Re  =  u  x.  /v  =  2.08  x  105 

A,  °°  K 


and  the  roughness  Reynolds  number  is 


Rek  =  550, 


defined  in  terms  of  the  protuberance  height  and  the  local  undisturbed  velocity 
at  x  =  xk,  z  =  k. 

Tani  took  hot  wire  surveys  of  the  protuberance  flow  field  at  various  x  to 
study  the  manner  in  which  the  streamwise  vortices  redistribute  momentum  in  the 
boundary  layer  and  affect  transition.  His  measurements  taken  at  (x-xk)/k  =  10 

are  used  in  the  present  study  to  construct  a  plane  of  data  for  "initializing" 
the  marching  schemes  used  to  solve  Equations  (14)- (16) .  From  the  limited  data 
taken  by  Tani,  distributions  of  velocity  components  which  fit  these  data  must 
be  constructed  in  order  to  generate  the  initial  plane  flow.  This  is  accom¬ 
plished  by  combining  the  available  data  with  theoretical  considerations,  so 
that  U,  V,  and  W  are  prescribed  everywhere  at  (x-xk)/k  =  10.  This  process  is 
described  in  the  next  section. 

In  addition  to  velocity  fields.  Reference  3  presents  data  for  measure¬ 
ments  of  the  flow  perturbation  caused  by  the  protuberance  and  the  spacing  of 
the  vortices  for  downstream  distance  (x-xk)/k  =  244;  these  data  are  used  to 

test  the  calculations.  The  transition  process,  which  starts  at  x  =  xt,  begins 
further  downstream;  from  Reference  3,  (xt-xk)/k  =  4500.  z 


IV.  INITIAL  PLANE  OF  DATA 

The  experimental  data3  upon  which  the  velocities  in  the  initial  plane,  xj 

=  xk  +  10k,  are  based  are  shown  in  Figure  2.  The  u  and  v  distributions  are 

consistent  with  two  vortices,  symmetrical  about  y  =  0,  but  are  insufficient  to 
characterize  the  entire  flow  field.  They  must  be  supplemented  with  theoreti¬ 
cal  considerations  and  assumptions.  The  deviation  of  the  fitted  values  from 
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the  actual  velocity  field  might  be  large  for  sorre  y  and  z  because  of  the 
paucity  of  the  data. 

Figures  2 a  and  2b  furnish  data  at  only  three  values  of  z,  whereas  values 
of  Uj-Ug  are  given  in  Figure  2c  for  many  z  values.  The  curves  in  these  fig¬ 
ures  are  from  Reference  3,  representing  the  authors'  smoothing  of  the  data. 
Figure  2c,  together  with  Figure  6  of  Reference  3,  shows  that  the  shape  of  the 
U(Y)  curves  is  the  same  for  all  Z;  only  the  magnitude  changes.  Three  quanti¬ 
ties  are  defined  graphically  in  Figure  2a:  U1  and  U2  are  the  maximum  and 

minimum  values,  respectively,  of  U  off  the  plane  of  symmetry  at  fixed  values 

of  l\  s/2  is  the  maximum  value  of  Y  at  which  U  =  j  (Uj  +  U2)  and  is  thus  a 

measure  of  the  spacing  of  the  two  legs  of  the  horseshoe  vortex.  A  measure  of 
the  perturbation  due  to  the  horseshoe  vortex  is  -  U2;  for  Figure  2  it  is 

approximately  20%  of  the  local  Blasius  value  of  U. 

The  fits  for  U  and  V  in  the  initial  plane  are: 


Uj  -  UB  +  1.961  [1.519  n  fQ"  (1.519  n)3  [UT(Y)  -  Uy  (2.5)] 


where 


Vj  =  -1.217  n  Uj  F2  (1.558  n)  tan  3y  (Y), 

n  -  2/(2  Xj)1/2  =  [u=o/(2v  xx  )1/2]z, 

UB  =  fo'  M 


(17) 

(18) 

(19) 

(20) 


are  the  Blasius  variables.  The  process  which  led  to  these  fitting  functions 
is  described  below. 

The  functions  fQ'  (n),  n  f  "  (n),  and  F^  (n)  are  tabulated  in  Appendix 
A.  These  are  obtained  by  solving 


f  (n)  +  f  f  "  «  0 
0  0  0 

F  ''  +  F9‘  +  3.3873  f  1  F0  =  0, 

Z  0  C  0  2 

where  fQ  (0)  =  fQ'  (0)  =  0;  fQ'  (»)  =  l 
F2  (0)  =  F2  (»)  =  0. 


(21) 
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F2  is  an  eigenfunction  occurring  in  the  boundary-layer  perturbation  analysis 
of  Fox  and  Libby5. 

The  functions  Uy(Y)  and  3y(Y)  are  obtained  by  dividing  the  Y-axis  into 

segments  and  empirically  fitting  the  data  of  Figure  2  in  each  segment  with  a 
3rd-order  polynomial,  with  continuity  of  the  functions  and  first  derivatives 
required  at  tie-points. 

uj(y)  =  a0  +  aiY  +  a2y2  +  a3y3 

o  ,  (22) 

3y(Y)  -  bQ  +  b^Y  +  b2Y^  +  bgY^  (in  degrees), 

where  3y  is  fitted  at  z  =  0.18  cm.  The  imposed  symmetry  requires  that  we  need 

to  consider  only  Y  >  0.  Five  segments  were  used  here.  The  ai  and  b*  coeffi- 

cients  and  tie-points  are  provided  in  Appendix  B.  Regrettably,  the  boundary- 
region  approximation  was  not  under  consideration  when  the  cubic  spline  fits 
were  made  for  Uy  and  3yj  consequently  Uyy  and  Vyy  in  Eqs.  (14)  and  (15)  are 

discontinuous  at  the  tie-points. 

The  following  process  was  used  to  arrive  at  Eqs.  (17)  and  (18).  We 

imposed  the  plausible  assumption  that  as  either  Y  »  or  Z  ®,  U  approaches 
Blasius  flow  and  V  approaches  zero.  In  Eq.  (17)  the  Blasius  U  is  reached  at  Y 
=  2.5.  The  function  n  f  ' '  (n)  has  the  same  shape  and  limits  as  the  Uj-l^ 

data  in  Figure  2c,  which  measure  the  amplitude  of  perturbation  to  Blasius 
flow.  The  two  constants  in  Eq.  (17)  were  obtained  by  imposing  two 
conditions:  (i)  U  should  agree  with  Tani's  measurement  at  y  =  0.30  cm  and  z  = 
0.25  cm  (Figure  2a),  and  (ii)  the  peak  of  the  n  f  1 1  (n)  curve  should  occur  at 

z  =  0.24  cm,  the  peak  of  the  data  in  Figure  2c.  As  for  Eq.  (18),  Figure  2b 

indicates  a  reversal  in  phase  of  the  3  versus  y  profiles  between  z  =  0.18  cm 

and  z  =  0.33  cm;  the  function  F£(o)  provides  this  feature.  The  two  constants 

in  Eq.  (18)  were  obtained  by  forcing  the  analytical  expression  for  Vj  to  agree 

with  the  data  at  y  =  0.30  cm,  z  =  0.18  cm,  and  y  =  0.30  cm,  z  =  0.33  cm. 

Results  of  the  fitting  are  compared  with  Tani's  measurements  in  Figure  3, 
which  shows  spanwise  distributions,  and  Figure  4,  which  shows  variation  normal 
to  the  plate.  In  Figure  3b  only  the  data  for  y  >  0  and  z  =  0.18  cm  were  used 
to  determine  the  constants,  and  the  deviation  between  these  data  and  the 
fitted  curve  is  relatively  small.  Although  the  deviation  for  y  <  0  is  rela¬ 
tively  large,  it  is  apparent  that  the  scatter  in  the  data  is  also  large,  e.g., 
the  data  for  3  is  not  an  odd  function  of  y.  In  fitting  the  data,  however,  we 
required  the  symmetry  conditions  to  be  satisfied. 


5.  Fox ,  H.,  and  Libby,  P.A.,  "Some  Perturbation  Solutions  in  Laminar 

Boundary  Layer  Theory .  Part  2.  The  Energy  Equation",  J.  Fluid  Mech., 
Vol.  19,  1964,  pp.  433-451.  - 
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The  vertical  component,  W,  is  found  from  U  and  V  by  subtracting  Eq.  (14) 
from  [U  times  Eq.  (16)],  then  separating  out  the  resulting  3(W/U)/9Z  term  and 
integrating: 

Z 

W  =  U  Iq  UPx  "  UVY  +  VUY  "  Re  1  Uyy  '  UZZ^/lj2^  dZ*  (23) 

A  similar  expression  can  be  derived  using  Eqs.  (15)  and  (16).  For  uniform 
external  flow,  Px  =  0.  The  Y-  and  Z-  derivatives  are  evaluated  by  finite 

differences  at  the  mesh  points  in  the  Y ,Z  plane  w'ch  3-point  formulas,  and  W 
is  obtained  by  numerical  integration  of  the  right-hand  side  of  Eq.  (23).  The 
integrand,  J,  of  Eq.  (23)  is  indeterminate  but  finite  at  Z  =  0,  so  that  at  the 
first  level  of  mesh  points  off  the  plate, 

W  (Z=AZ)  =  i  AZ  U(AZ)  J ( AZ/2) . 

The  computation  fails  if  U  =  0  off  the  plate;  however,  this  has  never  occur¬ 
red.  The  boundary- layer  and  boundary-region  approximations  yield  different 

distributions  because  of  the  absence  or  presence  of  the  Re"1  UYY  term. 

In  Figure  5  the  velocity  vectors  are  shown  projected  into  the  initial 
plane  (as  viewed  looking  downstream  from  the  protjberance) .  The  U-  and  V- 
components  are  identical  in  both  approximations,  out  the  W-  components  are 
different.  A  mesh  point  lies  at  the  tail  of  each  velocity  vector.  The  bound- 
ary-layer  approximation  leads  to  two  counter- rotatirg  vortices  on  each  side  of 
the  plane  of  symmetry.  In  the  boundary-region  approximation,  there  is  one 
large  vortex  and  possibly  several  smaller  vortices  that  are  not  clearly  de¬ 
fined.  One  noticeable  difference  occurs  at  the  plane  of  symmetry,  Y  =  0;  in 
the  boundary-region  case  the  flow  is  toward  the  plate,  while  in  the  boundary- 
layer  case  the  flow  is  away  from  the  plate. 

The  differences  between  the  two  flow  patterns  in  Figure  5  suggest  that 
UYY,  and  therefore,  VyY  are  not  negligible  for  this  vortical  flow.  Figure  6 

compares  uyy  with  uzz  at  Z  =  1.0.  (Since  these  two  terms  are  dimensional, 

they  may  be  compared  directly,  no  stretching  factor  being  involved.)  This 
figure  shows  that  Uyy  is  8  times  as  large  as  uzz  in  some  parts  of  the  initial 

plane.  Thus  the  boundary- layer  approximation  is  violated  in  this  plane 
because  it  requires  uyy  «  uzz. 

That  the  boundary-layer  approximation  may  be  volated  can  be  anticipated 
from  the  large  gradients,  e.g. ,  Uy  shown  in  the  data  in  Figure  2a.  The  ques¬ 
tion  arises  whether  or  not  a  boundary  layer  calculation  should  be  attempted. 
It  was  decided  to  do  this  on  the  possibility  that  the  large  gradients  would  be 
smoothed  out  and  that  the  downstream  flow  could  be  calculated. 


V.  NUMERICAL  PROCEDURE  AND  BOUNDARY  CONDITIONS 

The  equations  were  solved  using  the  predictor-corrector  multiple-itera¬ 
tion  ( PCM I )  method,  which  has  two  advantages.  It  can  be  used  for  both  sets 
of  equations,  and  the  truncation  error  is  uniformly  second  order.  One 
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alternative  procedure  for  the  boundary  layer  equations  is  discussed 
later. 

The  numerical  solution  to  Eqs.  (14)-(16)  is  obtained  by  marching  down¬ 
stream  (X-di rection)  from  the  initial  data  plane.  Symmetry  boundary  condi¬ 
tions 


3U/3Y  =  V  =  3W/3Y  =  0  (24) 

are  applied  at  Y  =  0  and  Blasius  flow  is  prescribed  at  Y  =  2.5.  The  free- 
stream  velocity  components  (U  =  1  and  V  =  0)  are  prescribed  at  Z  =  100.  The 
numerical  procedure  is  based  on  a  predictor-corrector  multiple-iteration 
(PCMI)  technique  developed  by  Rubin  and  Lin6.  The  procedure  is  implicit  in 
the  Z-di rection.  All  flow  gradients  in  the  Y-di rection  are  approximated  by 
prediction  and  subsequent  correction  using  an  iterative  approach,  which  gives 
a  more  accurate  simulation  of  the  nonlinear  coupling  between  equations.  The 
solution  is  obtained  for  each  grid  column  by  inverting  a  tridiagonal  matrix. 

Appendix  C  gives  the  set  of  difference  equations  used  in  this  study.  The 
truncation  errors  for  interior  points  are  of  0  (Ax2,  AY2,  AZ2).  The  PCMI 
procedure  solves  the  momentum  difference  equations  for  Um+'*'  and  Vm+1,  using 

the  m-iterate  values  to  form  the  coefficients  of  the  nonlinear  terms  and 
approximate  Y-deri vati ves ,  where  m  denotes  the  iteration  level.  The  continu¬ 
ity  equation  is  then  solved  for  Wm+1. 


For  the  predictor  step,  or  first  iteration,  terms  in  Eqs.  (C-l)-(C-lO) 


with  superscript  zero  are  approximated  by  a  Taylor  series  to  0  (aX^): 


F° 


1+1,  .i,k  Fi,j,k  +  (Fi,j,k  "  Fi-l,j,k)  (Xi+l_Xi ) /(Xi  “Xi_  i  )> 


(25) 


During  the  first  X-step  an  extrapolation  of  0  (AX)  is  used.  After  extrapolat¬ 
ing  guesses  for  U°,  V°  and  W°  at  X -.i ,  Eqs.  (C-l)  and  (C-6)  are  solved  to  give 

U  and  V  .  The  calculations  start  at  the  plane  of  symmetry  and  work  outward 
in  Y.  The  (m+1 )-iterate  values  are  used  to  approximate  derivatives  in  the  Y- 
di  rection  as  soon  as  they  become  available;  see  Eqs.  (C-5)  and  (C-10).  Next 

Eq.  (C-ll)  is  solved  for  W1,  starting  at  Y  =  Z  =  0  and  sweeping  in  Z  and  then 
Y.  This  completes  the  first  iteration  cycle.  Subsequent  iterations  allow  the 
nonlinear  terms  to  be  corrected  and  the  symmetry  plane  boundary  conditions  to 
be  satisfied. 


6.  Rubin,  S.G.,  and  Lin,  T.C.,  "A  Numerical  Method  for  Three-Dimensional 
Viscous  Flow:  Application  to  the  Hypersonic  Leading  Edge",  J.  Comp. 
Fhys.,  Vol.  9,  1972,  pp.  339-364.  - 
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Numerical  studies  were  performed  to  determine  the  number  of  iterations 
required  to  give  satisfactory  convergence  of  the  solution.  These  indicate 
that  three  iterations  assure  convergence  of  U,  V,  and  W  to  at  least  four 
decimal  places,  which  we  consider  to  be  sufficiently  accurate.  Consequently 
at  least  three  iterations  were  used  for  subsequent  calculations.  The  number  of 
iterations  used  varied  with  AX,  which  ranged  from  0.010  near  the  initial  data 
plane  to  0.250  at  the  last  downstream  station.  The  value  of  AX  was  adjusted 
during  the  calculations  to  satisfy  the  linear  stability  criterion 

AX  <  AY |U/V | .  (26) 


Eq.  (26)  is  based  on  a  linear  stablity  analysis  for  interior  points,  including 
the  effect  of  multiple  iteration6.  A  41  x  41  (Y  x  z]  grid  was  used  for  each  X 
=  constant  plane  with  AY  =  0.025  for  0  <  Y  <  0.5,  AY  =  0.100  for  0.5  <  Y  <  2.5 
and  AZ  =  2.50  for  0  <  Z  <  100.  The  variable  spacing  in  Y  was  adopted  to  help 
resolve  the  large  crossflow  gradients  indicated  in  Figure  2  for  0  <  Y  <  0.5. 

Special  second-order  accurate  difference  approximations  for  3F/3 Y  and  3^F/3Y^ 
are  employed  at  Y  =  0.5  due  to  the  unequal  mesh  spacing,  where  F  represents  U 
or  V. 


VI.  RESULTS  AND  COMPARISON  WITH  EX3ERIMENT 

Computations  were  carried  out  for  the  vortical  flow  downstream  of  the 
protuberance  using  both  the  boundary- layer  and  boundary-region 

approximations.  The  boundary-layer  calculation  failed  when  it  predicted 
streamwise  flow  reversal;  the  boundary- region  calculation,  on  the  other  hand, 
was  carried  to  244  protuberance  heights  downstream  (Tani's  last  station) 
without  encountering  computational  difficulties. 

Crossflow  velocity  vector  plots  at  (x  -  xk)/<  =  16  are  presented  in 

Figure  7,  which  shows  the  initial  stages  of  the  breakdown  of  the  boundary- 
layer  calculation.  Large  gradients,  Wy,  are  predicted  at  Y  =  0.025  and  0.3 

which  were  not  present  in  the  initial  plane.  As  -he  calculation  proceeds, 
they  are  not  smoothed  out,  due  to  lack  of  crossflow  ciffusion.  The  positive  W 
along  the  plane  of  symmetry,  present  on  the  initial  plane,  is  still  in  evi¬ 
dence.  At  x  =  xk  +  20  k  the  W-component  is  larger  and  the  U-component  re¬ 
verses  direction,  as  if  streamwise  separation  had  occurred.  Since  there  is  no 
pressure  gradient  in  this  flow,  the  usual  cause  of  separation,  an  adverse 
pressure  gradient,  cannot  be  invoked.  We  only  note  that  (i)  for  both  the 
initial  plane  and  for  (x-xk)/k  =  16,  W  >  0  in  the  neighborhood  of  Y  =  0,  and 

(ii)  Uz  >  0  up  to  the  point  of  flow  reversal;  thus  the  term  WUZ  could  play  the 

role  of  an  adverse  pressure  gradient.  At  any  rate  the  boundary-layer  calcula¬ 
tion  cannot  proceed  beyond  x  =  x.,  +  20  k. 

For  the  boundary-region  calculation  the  crossflow  velocity  vector  plot  in 
Figure  7  shows  that  some  anomalies  that  existed  in  che  initial  plane,  caused 
by  discontinuities  in  second  derivatives  with  respect  to  Y,  have  been  smoothed 
out.  The  presence  of  the  main  vortex  is  clearly  seen.  In  contrast  to  the 
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boundary  layer  results,  W  <  0  near  Y  =  0  at  this  plane  as  it  was  in  the  ini¬ 
tial  plane,  and  the  gradients  are  smaller,  the  largest  Wy  for  the  boundary 

region  case  being  20%  of  that  for  the  boundary  layer  case.  These  differences 
between  the  results  of  the  two  approximations  are  felt  to  be  related  to  the 
success  of  the  boundary  region  calculation. 

Figure  8  shows  the  downstream  variation  of  the  maximum  value  of  the 
streamwise  vorticity,  c  (=  9W/9Y  -  9V/9Z),  evaluated  at  the  wall.  This  is  a 

measure  of  the  strength  of  the  horseshoe  vortex.  Results  from  the  boundary- 
layer  calculation  are  terminated  where  it  breaks  down.  The  streamwise  vortic¬ 
ity  from  the  boundary-region  calculation  is  essentially  zero  at  about  150 
protuberance  heights  downstream. 

With  regard  to  crossflow  velocity,  the  boundary-region  calculation  pre¬ 
dicts  that  the  flow  inclination  at  (x  -  xk)/k  =  244  is  less  than  0.2°.  Tani 

could  not  detect  any  crossflow  velocity  component  at  244  protuberance  heights 
downstream.  In  fact,  0.2°  is  within  scatter  of  his  data. 

Figure  9  shows  the  streamwise  variation  of  the  horseshoe  vortex  spacing, 
s.  The  variation  predicted  by  the  boundary-region  approximation  is  seen  to 
agree  qualitatively  with  Tarn's  measurements.  At  the  last  downstream  station 
the  predicted  spacing  is  about  25%  larger  than  the  measurement.  The  trend  of 
the  calculated  s  is  the  same  as  that  of  the  data. 

Figure  10  shows  the  variation  of  the  maximum  velocity  difference,  (Uj- 
Uo)  ,  .  The  quantity  Ui-Uo  is  evaluated  in  each  crossflow  plane  and  the 
largest  value,  (U j-U2)max  is  plotted.  The  boundary-region  calculation  pre¬ 
dicts  that  the  perturbation  in  U  grows  to  ~  0.19  at  x  -  x^  +  25  k  and  then 

begins  to  decay.  The  calculated  perturbation  decays  faster  than  the  measured 
value  up  to  (x  -  xk)/k  *  200;  thereafter  the  calculated  value  is  essentially 

constant  and  is  38%  of  the  experimental  value  at  244  protuberance  heights 
downstream.  Note  that  the  perturbation  in  U  is  still  present  at  (x  -  xk)/k  = 

244,  although  the  crossflow  velocity  has  decayed  to  nearly  zero  and  the  maxi¬ 
mum  streamwise  vorticity  at  the  wall  is  already  zero.  These  conclusions  apply 
also  to  the  experimental  data. 

The  failure  of  the  boundary  layer  calculation  would  seem  to  be  related  to 
the  lack  of  a  mechanism  for  smoothing  out  or  diffusing  the  large  crossflow 
gradients  in  the  initial  data.  However,  care  must  be  taken  in  drawing  conclu¬ 
sions  from  the  results  of  these  numerical  solutions  because  the  failure  of  the 
PCMI  boundary- layer  calculation  may  be  a  product  of  the  numerical  scheme 
and/or  the  inadequacy  of  the  boundary  layer  model  itself.  An  indication  might 
be  obtained  by  trying  other  schemes.  Rather  than  try  schemes  similar  to  the 
PCMI,  which  would  probably  fail  in  the  same  way,  we  used  a  scheme  that  adds 
support  to  the  statement  made  in  the  first  sentence  of  this  paragraph.  The 
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scheme  chosen  was  the  one  in  Reference  7,  called  Scheme  D,  on  the  basis  of 
zone-of-dependence  sensitivity.  For  details  see  Reference  7.  The 


truncation  errors  for  this 

with  the  PCMI  truncation  error, 

2 

the  additional  error  AY  /AX.  The  latter 
introduced  by  the  X-difference  equation. 


scheme  are  0(AX,  AY^ ,  AZ^,  AY2/AX).  In  contrast 

this  scheme  is  only  "irst  order  in  AX  and  has 

from  artificial  diffusion 
standard  truncation  error 
two  artificial  viscosity 

introduced  into  the  finite- 


analysis  techniques, 

terms,  a  Uyy  and  a  Vyy , 

difference  equations8, 
form  as  the  crossflow 
These  terms  have  the 
gradients, 
layer  PCMI 


it  can  be  shown 


arises 

Using 

that 


where  a  =  AY^/AX,  are 


These.  two  artificial  viscosity  terms  have  the  same 
diffusion  terms  in  the  boundary- region  approximation, 
effect  of  smoothing  out  the  initial  large  crossflow 
and  the  calculation  proceeds  beyond  the  point  where  the  boundary 
calculation  failed.  Apparently  the  Scheme  D  boundary-layer 


difference  equations  behave  in  a  manner  simila'r  to  the  boundary- region  approx 
imation.  By  increasing  AX  from  0.01  at  the  initial  plane  to  0.50  far  down¬ 
stream,  the  calculation  proceeded  to  X  =  244.  However,  the  results  compared 
poorly  with  the  data.  Figure  11  shows  the  values  of  s  from  this  scheme  to- 
gether  with  the  boundary-region  result.  The  divergence  between  the  results, 
as  X  increases,  may  be  caused  by  the  fact  that  a  is  50  times  smaller  down¬ 
stream  than  near  the  initial  plane,  and  hence  the  artificial  diffusion  is 
small,  or  perhaps  because  AX  is  50  times  larger  thar  near  the  initial  plane 
and  thus  produces  unacceptable  truncation  errors.  The  conclusion  is  that  this 

boundary- layer  calculation  fails  but  in  a  different  sense  from  the  PCMI  fail- 
ure. 


VII.  CONCLUSIONS 

In  this  study  we  have  shown  that  the  three-dimensional  boundary-region 
approximation  can  be  successfully  used  to  describe  the  flow  downstream  of  a 
cylindrical  protuberance  immersed  in  a  laminar  flat  plate  boundary  layer.  The 
persistence  of  the  streamwise  vorticity  downstream  of  the  obstacle  is  predict¬ 
ed  in  agreement  with  experimental  observations.  We  find  that  the  streamwise 
vorticity  persists  more  than  100  protuberance  heights  downstream  and  that  the 
initial  perturbations  in  U  decay  more  slowly  than  the  vorticity.  Significant 
perturbations  in  U  are  still  present  244  protuberance  heights  downstream  of 
the  obstacle  in  the  boundary-region  calculation.  This  is  quite  remarkable 


7 *  Kitchene  Jr.,  C.W.,  Sedney,  R.,  and  Gerber,  N.,  "The  Role  of  the  Zone  of 
Dependence  Concept  in  Three-Dimensional  Boundary-Layer  Calculations" ,  U.S. 
Army  Ballistic  Research  Laboratory/ ARRADCOM  Report  No.  1821,  Aberdeen 
Proving  Ground,  MD,  August  1975.  AD  A016896. 

3‘  "0n  Artificial  Viscosity",  J.  Comp.  Phys.,  Vol.  10,  October 

1972,  pp.  169-184.  - 
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considering  that  we  have  marched  more  than  a  hundred  boundary- layer  thick¬ 

nesses  downstream  to  reach  this  position.  The  experimental  data  show  the  same 
qualitative  trend. 

The  three-dimensional  boundary-layer  approximation,  on  the  other  hand,  is 
not  successful  in  describing  the  flow  field  downstream  of  the  protuberance. 
It  fails  due  to  the  large  crossflow  gradients  present  in  the  protuberance  flow 
field.  These  gradients  lead  to  the  reversal  of  the  streamwise  velocity  compo¬ 
nent,  U,  during  the  marching  process. 

It  is  expected  that  the  quantitative  agreement  between  the  boundary- 

region  calculations  and  Tani's  experimental  data  would  be  improved  if  the 

initial  data  were  fitted  using  4th-order  (or  higher)  polynomials  so  that  UYY 

and  VYY  were  continuous.  Accuracy  could  be  further  improved  by  using  smaller 

grid  sizes  during  the  marching  process.  Unfortunately,  it  is  not  possible  to 
quantify  the  expected  improvement.  These  refinements  may  not  be  warranted  in 
the  present  problem,  considering  the  sparse  experimental  data  available  for 
constructing  initial  data. 

Whether  or  not  the  boundary- region  approximation  employed  here  is  the 

most  adequate  model  for  this  flow  field  cannot  be  established  unequi vocably  by 
this  investigation.  In  our  version  we  assumed  that  the  normal  pressure  gradi¬ 
ent,  3P/3Z,  is  zero.  It  is  possible  that  dropping  this  assumption  would 
improve  the  quantitative  agreement  with  experiment. 
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Figure  1.  Sketch  of  Protuberance  in  Boundary  Layer  in  Tani 
Experiments  =  0.39,  Re  =  2.08  x  105,  Re.  = 
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and  (b)  Flow  Inclination  in  the  x,y  Plane  at  Two  Heights;  (c)  Experimental 
Distribution  Across  Boundary  Layer  of  Peak-Valley  Velocity  Difference  (U 


Figure  3.  Analytical  Spanwise  Distributions  of  (a)  u/ua  and  (b) 

3  =  tan-1  (v/u),  Eqs.  (17)  and  (13),  at  Fixed  Heights 
in  the  Initial  Data  Plane,  and  Comparison  with  Tani  Data. 
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Figure  4.  Analytical  Distributions  of  (a)  u/u^  and  (b)  U]  -  U2  Through  the  Boundary  Layer  in 
the  Initial  Data  Plane,  and  Comparison  with  Tani  Data. 
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Figure  5.  Crossflow  Plane  (y,z)  Velocity  Vector  Projections  in  Initial 
Data  Plane  at  x  =  x.  +  10k 


y  (cm) 


Figure  6.  Spanwise  Distribution  of  Second-Derivatives  of  Streamwise 
Velocity  Component  in  Initial  Data  Plane 


26 


o 

II 


g 

u 


o 


QC 

LU 

$ 


> 

QC 

2 

Z 

o 

CD 


£ 


VO 

+ 

-X 

X 

II 


27 
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Figure  7.  Crossflow  Plane  (y,z)  Velocity  Vector  Projections  at  x 
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Figure  8.  Streamwise  Variatior  of  Maximum  Vorticity  at  Wall 
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Figure  9.  Streamwise  Variation  of  Horseshoe  Vortex  Spacing 
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Figure  10.  Streamwise  Variation  of  Maximum  Peak-Valley  Velocity 
Difference  Induced  by  Horseshoe  Vortex 


Figure  11.  Comparison  of  Horseshoe  Vortex  Spacing:  3-D  Boundary 
Layer  (Scheme  D)  vs  Boundary  Region 
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LIST  OF  SYMBOLS 


P 

P 

Re 


Rek 


s 

u  ,v,w 


Uo> 

U,V,W 
Ui  9^2 
UB 

UI»VI 

UT 


x»y,z 


xk 


X  ,Y  ,Z 
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functions  employed  in  formulation  of  initial  plane  data,  Eq.  (21) 
height  of  cylindrical  protuberance  (=0.25  cm  in  expt.) 
reference  length  for  nondimensional izinq  (=  1  cm) 

O 

pressure  [g/(cm  s  )] 

=  P/(PUoo2) 

=  pu^L/p  =  Reynolds  number  (=  3467  in  expt.) 

=  puj</y  =  roughness  Reynolds  number  (=  550  in  expt.) 

=  pu^x^/y  (=  2.08  x  105  in  expt.) 

2  times  value  of  y  at  which  U  =  j  (Uj  +  U2),  Figure  2a  [cm] 

velocity  components  in  the  directions  of  x,  y,  and  z,  respectively 
[cm/s] 

free  stream  velocity  (=  520  cm/s  in  expt.) 

1  /2 

=  u/u^,  v/u^,  Re  w/u^,  respectively 

maximum  and  minimum  values  of  U,  respectively,  (Figure  2a) 

Blasius  flat  plate  boundary  layer  velocity,  Eq.  (20) 

U  and  V  at  X  =  Xj ,  plane  of  initial  data 

U-velocity  function  obtained  by  fitting  Tani  data,  Eqs.  (17)  and 

rectangular  coordinates:  streamwise,  spanwise,  and  normal  to  wall, 
respectively  (Figure  1)  [cm],  origin  at  leading  edge  of  flat  plate 

x-coordinate  of  cylindrical  protuberance.  Figure  1  (=  60  cm  in 
expt. ) 

=  x/L,  y/L,  Re1/2  z/L,  respecti vely 

X-coordinate  of  initial  data  plane  (=  62.5  in  expt.) 

function  describing  spanwise  VT/UT,  obtained  by  fitting  Tani  data, 
Eqs.  (18)  and  (22)  11 

boundary  layer  thickness.  Figure  1  [cm] 
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LIST  OF  SYMBOLS 
(Continued) 


AX ,aY  ,AZ  increments  in  X,  Y,  and  Z  in  finite-difference  equations,  Appendix 


=  0  for  boundary  layer,  =  1  for  boundary  region,  Eqs.  (14),  (15), 
and  (23) 


5x 

n 

u 

v 

P 


streamwise  vorticity  [ 3W/9Y  -  3V/3Z) 

=  Z/(2  X,)V2,  Eq.  (19) 

viscosity  of  air  [g/(cm  s)] 

kinematic  viscosity  o~  air  (=  0.15  cm^/s  in  expt) 
density  of  air  [g/cm^] 


Subscript 


I 

k 

T 


indices  identifying  grid  points,  Eq.  (25)  and  Appendix  C 

initial  plane  of  data 

front  of  cylindrical  protuberance 

Tani  measurement 

free  stream 


Superscri pt 

m  designator  of  iteration  level 
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APPENDIX  A:  TABLES  OF  fQ' ,  n  fQ"  ,  and  F2 


n 

V 

_  £  1  1 
n  fo 

F2 

0.00 

.0000 

.0000 

.0000 

0.25 

.1174 

.1173 

.2494 

0.50 

.2342 

.2325 

.4906 

0.75 

.3493 

.3408 

.7032 

1.00 

.4606 

.4344 

.8574 

1.25 

.5656 

.5044 

.9223 

1.50 

.6615 

.5427 

.8776 

1.60 

.6967 

.5480 

.8286 

1.75 

.7458 

.5448 

.7254 

2.00 

.8167 

.5113 

.4948 

2.25 

.8736 

.4490 

.2354 

2.50 

.9168 

.3687 

.0013 

2.75 

.9479 

.2831 

-.1671 

3.00 

.9691 

.2031 

-.2543 

3.25 

.9826 

.1363 

-.2691 

3.50 

.9907 

.0855 

-.2353 

3.75 

.9953 

.0501 

-.1795 

4.00 

.9978 

.0275 

-.1226 

4.25 

.9990 

.0141 

-.0759 

4.50 

.9996 

.0068 

-.0430 

4.75 

.9998 

.0031 

-.0224 

5.00 

.9999 

.0013 

-.0108 

5.25 

1.0000 

.0005 

-.0048 
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APPENDIX  B:  POLYNOMIAL  COEFFICIENTS  IN  EQ.  (22) 
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APPENDIX  C:  PCMI  FINITE  -  DIFFERENCE  EQUATIONS 


Uncoupled  linear  finite-difference  equations  are  constructed  with  the 
PCMI  method  by  approximating  U,  V  and  W  in  Eqs.  (14)  and  (15)  as  known 
constants,  specified  from  either  the  extrapolated  value  or  the  value 
calculated  at  the  previous  iteration  level.  The  coupling  between  equations 
and  the  nonlinearities  of  individual  terms  are  approximated  by  the  multiple- 
iteration  process  which  updates  U,  V  and  W  at  each  subsequent  iteration 
level.  The  X-momentum  difference  equation  is  given  by  the  tri diagonal  system 


m+1 

i+1 ,  j  ,k-l 


U 


ntf-1 

i+l,j,k 


U 


m+1 

i+l > j  >k+l 


(C-l) 


where 

aij  ■  -Ki,s, k  +  Wi,j,k)/(8AZ>  -  1/(2az2) 

blj  *  (UHl,j,k  +  1Ji,j,k)/(2AX)  *  iv2> 

clj  =  “aj j  -  1/(AZ  ) 


(C-2) 

(C— 3 ) 

(C-4) 


*U  “  u1.J.k  <u1*l,j,k  +  Ui,j,k)/(2AX) 

+  <ui>j.k+i-2Ui,j,ktui,j,k-i>/(2Az2) 

-  (“f.l.j.k  +  Mi'.j.k)(ui,j,w  -  U1.j,k-1>/<8AZ> 

_(Vm  +  V  1  (um  -  u"*1  +  u 
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with  e  =  0  in  the  3D  boundary  layer  approximation  and  e  =  1  in  the  3D  boundary 
region  approximation.  The  Y-momentum  difference  equation  is  given  by  the 
tridiagonal  system 


where 


,  x/nH-l  ,  .  ,,m+l  ,  .  ,,m+l  _  . 

a2j  Vi+l,j,k-l  b2j  ^i+1 ,  j ,k  C2j  Vi+l,j,k+l  “  d2j 


(C-6) 


a2j  =  al j  ’ 

(C-7) 

b2j.  =  blj# 

(C-8) 

c2.  =  clj# 

(C-9) 

d2j“Vi,j,k  «UT+l.J.k  +  U1,J.k>^2AX> 
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"  (Wi+l,j,k  +  Hi.J,k»V1J,k+l  “  Vi,j,k-l)/(8AZ) 

(c-io) 
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-2Vf,j,k  +  vU-l,k)/(2RelY2)- 
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The  continuity  difference  equation  is  given  by 


Clj>k  ■  “wd.k-1 + 
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fference  approximation  for  8U/3X  is  based  on  a 

second- 

order  accurate  backward-difference  formula  for  unevenly-spaced  points 
involving  three  X-levels;  the  first  step  in  X  off  the  initial  data  plane 
requires  a  difference  formula  involving  only  two  X-levels. 
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